Gravitational-wave bursts and stochastic background from 
superfluid vortex avalanches during pulsar glitches 
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ABSTRACT 

The current-quadrupole gravitational-wave signal emitted during the spin-up phase 
of a pulsar glitch is calculated from first principles by modeling the vortex dynamics 
observed in recent Gross-Pitaevskii simulations of pinned, decelerating quantum con- 
densates. Homogeneous and inhomogeneous unpinning geometries, representing creep- 
and avalanche-like glitches, provide lower and upper bounds on the gravitational wave 
signal strength respectively. The signal arising from homogeneous glitches is found 
to scale with the square root of glitch size, whereas the signal from inhomogeneous 
glitches scales proportional to glitch size. The signal is also computed as a function of 
vortex travel distance and stellar angular velocity. Convenient amplitude scalings are 
derived as functions of these parameters. For the typical astrophysical situation, where 
the glitch duration (in units of the spin period) is large compared to the vortex travel 
distance (in units of the stellar radius), an individual glitch from an object 1 kpc from 
Earth generates a wave strain of 10-24[(Aa;/w)/10^'^](a;/102 rads~^)3(Ar/10"2 ^^)-i^ 

where Ar is the average distance travelled by a vortex during a glitch, Auj/lj is the 
fractional glitch size, and w is the pulsar angular velocity. The non-detection of a signal 
from the 2006 Vela glitch in data from the fifth science run conducted by the Laser In- 
terferometer Gravitational- Wave Observatory implies that the glitch duration exceeds 
~ 10^* ms. This represents the first observational lower bound on glitch duration to 
be obtained. 
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1 INTRODUCTION 

Glitches are stochastic spin-up events that interrupt the 
smooth electromagnetic spin down of a pulsar. They 
are commonly attributed to the simultaneous, transi- 
tory unpinning and outward moti on of between 10^ and 
10^^ quantised s uperfl u id vortices IIAnderson fc Itoh 19751: 
Anderso n et al.l 11 982: 'Alpar fc Sauls' '1988: Cheng et al.1 
1988: Link fc Epstein 1996: Warszawski fc Melatos 200i : 
Melatos fc Warszawskil 120091 ) . whi ch pin metastably to nu- 
clei in the crystalline inner crust (jPines et al.lll98Cl ) except 
during a glitch. The internal superfluid reorganises its ve- 
locity field nonaxisymmetrically during such an event, driv- 
ing a time- varying current quadrupole moment, which emits 
gravitational radiation. 

Pulsar glitches are sources of both burst and continu- 
ous gravitational waves. Previous theoretical and observa- 
tional studies identified the post-glitch recovery phase, dur- 
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ing which the viscous fluid component comes into corotation 
with the crust, as a continuous- wave source, with particular 
attention paid to how the shear viscosity, charged fluid frac- 
tion and mutual friction parameter of neutron stars may be 
extracted from the p rospect i ve signal ^ Andcrsson & Come: 
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frequency-time maps of cross-power between two spatially 
separated terrestrial detectors, was proposed recently to 
identify 'long transient' ev ents lasting from sec onds to weeks 
like post-glitch relaxation (IThrane et al.ll201Cl) . along with a 
multi-detect or Bayesian algo rithm using J^-statistic ampli- 
tude priors (jPrix et al.ll201lh . 

The burst signal may be split into two categories. The 
flrst categor y includes radiatio n from inertial r-mode os- 
cillations ( |Glampedakis fc And ersson 2009) or quasi-radial 
acoustic and inertial mode s involving one or two fluid com- 
ponents (jSiderv et aLlbOlOT ). The strongest signal is achieved 
when the velocities of the neutron superfluid and the charged 
fluid oscillate in anti-phase. For a detector with a spec- 
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tral noise density of lO'^^Hz"^'^ at 103 Hz, the maximum 
signal-to-noise ratio for a glitch from a pu lsar at the distance 
of Vela is 0.9 (lAnde rsson fc Comeiil2001h . The second cate- 
gory of signal, resulting from sudden superfluid deceleration 
when many vortices unpin simultaneously and move radially 
outwards, is discussed in this paper. 

Many other theories have been proposed for how this 
happens. Most glitch mechanisms appeal to the storage of 
angular momentum in the superfluid interior, followed by 
its rapid release to the crust in disc rete events via the sud- 
den unpinni g of superf luid vortices (I Anderson fc Itohll 19751 : 
lAlpar et al.l ll984: Hask ell et al.l2011^ ■ Many ideas have been 
proposed for how this happens , the main point o f distin ction 
being the unpinning trigger. lAndersson et al.l (1200J) sug- 
gested that an instability between the viscous and inviscid 
fluid components, governed by the two flu ids' relative ve- 
locity, can trigger a glitch. iLink fc EpsteinI (|l996. ) invoked 
a sudden change in the coupling of the two fluid compo- 
nents, following the injection of heat into the interior, for 
example after a sudden rupture in the star's crust, known 
as a star quake. Star quakes, in which the crust fails and 
abruptly lowers its moment of inertia, are also invoked as a 
stand-alone glitch theor y for some objects (jRudermanll 19691 : 
iMiddleditch et al]|2006l '). although star quakes on their own 
cannot account f or the regular, large glitches observed in 
pulsars like Vela (jCrawford fc Demianskil l2003'l. 

Gravitational waves result from non-axisymmetric re- 
arrangements of mass and/or momentum field, so differ- 
ent glitch theories predict gravitational wave signals of 
different form and strength. Therefore it is not possible 
to predict the form of the glitch gravitational wave sig- 
nal independently of a glitch model, even though one can 
place li mits on its maximum str ength from energy argu- 
ments (jAndersson fc ComeI^l2001^ . In this paper, we con- 
sider specifically the self-organised, catastrophic unpinning 
model. The model is inspired by the statistics of glitches 
[scale-free (power-law) sizes and independent (Poissonian) 
waiting times], which point to an underlying process in- 
volving the collective behaviour of many individual elements 
l|Melatos et al.ll2008l ). 

A Bayesian search technique (jClark et al.l 120071 ) was 
recently applied to data collected during the fifth science 
run of the Laser Interferometer Gravitational- wave Observa- 
tory (LIGQj) to look for gravitational radiation from quasi- 
normal mode oscillations excited by a glitch in the Vela 
pulsar. No signal was detected. The glitch, detected elec- 
tromagnetically on 12 August, 2006, at the Hartebeesthoek 
Radio Observatory in South Africa, occurred during a pe- 
riod of five and a half hours during which the two LIGO de- 
tectors at Hanford returned high-quality, contiguous data. 
Data from the third detector at Livinginston were not in- 
cluded in the search due to seismic noise. The search re- 
turned an upper bound on the intrinsic strain amplitude of 
between 6.3 x lO"'^^ an d 1.4 x 10~'^°, depen di ng on the wave 
numb er of the mode (jAbadie et al.l l201ll ). iHavama et al.l 
(|2008l ) employed a Monte-Carlo method to study the de- 
tection efficiency of possible gravitational waves triggered 
by the 2006 Vela pulsar glitch. Glitches in Vela recur quasi- 
periodically, at intervals of approximately three years, so 
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there is a case for altering the observation schedule to en- 
sure that a world-wide detector network is operational for 
the next glitch [th e most recent Vela glitch occurred in Au- 
gust, 2006 iFlanagan fc Buchneij (fJOOQ )]. 

Radio telescope timing experiments on pulsar glitches 
reveal power-law-distributed sizes and exponential wait- 
ing times. These statistics point to an underlying collec- 
tive process that transitions b etween metastable s tates via 
vortex unpinning avalanches (jMelatos et al.M2008r ). Quan- 
tum mechanical simulations confirm local and global knock- 
on mechani sms for triggering and sustaining unpinning 
avalanches IjWarszawski et al.ll2012h . In many models, the 
avalanches are spatially localised, so the glitch process is 
intrinsically nonaxisymmetric. 

In this paper, we calculate the gravitational wave signal 
from vortex avalanches during the spin-up phase of a pulsar 
glitch. Section sec:GPE summarises the first-principles evi- 
dence from quantum mechanical simulations, that a broad- 
band gravitational wave signal arises from non-axisymmetric 
superfiuid vortex rearrangement. The simulations consist of 
numerical solutions of the non-linear Schrodinger equation, 
known as the Gross-Pitaevskii equation, which describes the 
dynamics of a Bose-Einstein condensate. Such simulations 
are limited to small systems containing ~ 10^ vortices. In 
Sec. 13] onwards we incorporate the generic behaviour ob- 
served in the Gross-Pitaevskii simulations into an approx- 
imate analytic theory of vortex motion, which is valid for 
realistically large numbers of vortices. Sections [3] and 2] eval- 
uate the current quadrupole moment generated by spas- 
modic vortex motion in an idealised neutron star geometry. 
In Sec. [S] we calculate the gravitational wave signal from 
the motion of a single vortex. The signal arising from a vor- 
tex avalanche is calculated in Sec. |6l including considera- 
tion of the avalanche opening angle (Sec. l6.l"|) . vortex speed 
(Sec. 16. 4p and glitch size (|6.5|l . Criteria are given for detect- 
ing a gravitational wave burst from an individual glitch. In 
Sec. [3 we summarise our results and discuss critically the 
underlying assumptions of the calculation. For completeness, 
we also explain in the appendix how the theoretical ma- 
chinery in Sec. O-El can be applied to calculate rigorously 
the glitch contribution to the stochastic gravitational-wave 
background. The signal-to-noise ratio is small for current- 
and next-generation detectors, if only Galactic pulsars are 
considered, but it may be significant when extragalactic ob- 
jects are included, an interesting avenue for future work. 



2 MOTIVATION: QUANTUM MECHANICAL 
SIMULATIONS 

Let us begin by confirming that, at the microscopic level, a 
quantum condensate (such as the neutron superfluid in the 
pulsar interior) in a decelerating container undergoes spas- 
modic vortex reorganisation, which emits a bursty broad- 
band gravitational wave signal. Such systems are tradition- 
ally studied by solving the Gross-Pitaevskii equation (GPE), 
a non-linear Schrodinger equation of the form 



(1) 



where i/> is the two-dimensional order parameter of a zero- 
temperature Bose-Einstein condensate, from which the con- 
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Figure 1. Top: Greyscale snapshots of superfluid density |-(/)p (wliere i/i is tiie two-dimensional complex order parameter) at t = 0, 
100, 200, 560 and 810 (in arbitrary units) from left to right respectively. The darker dots, and the lighter dots that are not part of 
the rectangular array, are vortices; lighter dots in the array are unoccupied pinning sites. Middle: Squared modulus of the gravitational 
wave strain in the cross polarisation |/ixP as ^ function of time t {black curve plotted against left vertical axis) and angular velocity 
u}(t) of the crust as a function of time {grey curve plotted against right vertical axis) for the simulations described in Sec. [2] (further 
detail is given in IWa rszawski & Mclatos 2011). Bottom: Frequency-time spectrogram of \hx P calculated in time windows of width 12.3. 
Simulation parameters: Nc /Ic = 10~^°, Vp = 16.6, rj = 1, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 0.0025, ujq = 0.8 (definitions of 
these quantities appear in IWarszawski fc Melatoal201M . 



densate density, \tp\'^, is derived, and |i/'P'(/' is the quan- 
tum pressure term which governs the structure of vortex 
cores. V is the potential, which comprises the container 
that confines the condensate, as well as a grid of spikes, 
which act as pinning sites for vortices, since the overlap of 
the empty vortex core with the region emptied of conden- 
sate by the high potential re sults in a low-energy configu- 
ration (|Avogadro et al.l 120081 1. V also contains a term pro- 
portional to the chemical potential. Initially, when the con- 
tainer is accelerated, vortices form and pin to the pinning 
sites. The initial pinning positions depend on the number 
of vortices and t h e number and strength of pinning sites 
JSato et al.l l2007l: iMink et al.l l2009l : iGoldbaum fc Muelled 
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( Jackson fc Barenghill20oi )^ 



201 ih and exhibit hysteresis 



Previous numerical studies of a vortex lattice, pinned to 
a decelerating container reveal that vortices hop spasmodi- 
cally betwee n pinning sites as they move towards the con- 
tainer edge (jWarszawski fc Melatod |2011^ . Computational 
expense limits the system size to tens of vortices, so the re- 
sults do not translate directly to a realistic pulsar. Instead, 
we use the qualitative behaviour observed as a basis for an- 
alytic extrapolation in later sections, when calculating the 
gravitational radiation from real glitches. 

Our simulations solve the Gross-Pitaevskii equation 
(GPE) for a zero-temperature condensate in a rotating con- 
tainer. For simplicity, we model an infinite vertical column of 
superfluid (strictly, its two-dimensional cross-section) in the 
presence of a rectangular grid of pinning sites. A summar y 
of the technical details is given in IWarszawski et al.l l|2012h . 
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Figure 2. Constant-t cross-sections of the spectrogram of |/ixP graphed in Fig. [T] at i = 551.9, 564.2, 576.5, 797.2, 809.5 and 821.8 
(black dotted, black solid, black dashed, grey dotted, grey solid and grey dashed curves respectively), showing spectral power as a function 
of frequency on log- linear axes (arbitrary units). 



The initial state of the system, depicted in the top left image 
of Fig. [T] (a greyscale snapshot of superfiuid density), com- 
prises 12 pinned vortices (darker dots) and 38 unoccupied 
pinning sites (lighter dots). An external spin-down torque is 
then applied to the container, imitating the electromagnetic 
spin-down torque acting on a pulsar. Pinning prevents vor- 
tices from moving radially outwards, fostering an angular 
velocity mismatch between the superfiuid and container. In 
response, vortices hop between pinning sites in spasmodic 
bursts, releasing discrete parcels of angular momentum to 
the container. The images in the top panel of Fig. [T] depict 
the superfiuid density at i = 0, 100, 200, 560 and 810 (the 
external torque is turned on at f = 0; time is quoted in arbi- 
trary simulation units). We can see that vortices move from 
their initial positions. For example, the interstitial vortex in 
the bottom right quadrant at i = annihilates against the 
container wall by t = 560. 

From movies of the time evolution of the superfiuid den- 
sity, we see that the typical vortex trajectory is not simple. 
As the vortex prepares to unpin, it migrates slowly to the 
outer edge of its pinning site. It then accelerates off the pin- 
ning site, before spiralling into the new pinning site if and 
when it repins. In later sections, we approximate the radial 
speed as a parabolic function of time to facilitate analytic 
calculation s. The motion of individu al vortices is described 
in detail in I Warszawski et al.l (J2012f ). 

Each time a vortex moves, the superfiuid velocity field 
adjusts in response. Each vortex generates a solenoidal ve- 
locity field, whose magnitude is inversely proportional to 
distance from the vortex core. The total velocity, Vg, at any 
point is the vector sum of the contributions from each vor- 
tex. Time-varying nonaxisymmetries in Vg generate a time- 
varying curre nt quadrupole moment, resultirig in gravita- 
tional waves (|Thornelll980l : lMelatos fc Perahall2"oTol ). 



We graph the amplitude squared of the wave strain 
in the cross polarisation, j/ixP (in arbitrary units), whose 
mathematical definition is formalised in Sec. O as a black 
curve in the middle panel of Fig. [l] The grey curve in 
the same plot tracks the angular velocity of the container, 
Lu. Note that all calculations are performed in the refer- 
ence frame that co-rotates with the container and pinning 
grid. Episodes of non-zero |/ixP are clearly accompanied 
by crustal spin-up events, i.e. glitches, when cj increases 
abruptly. Notably, the maximum \hxf during a glitch is 
not a monotonic function of glitch size (measured as Au/oj, 
the fractional, impulsive change in uj during a glitch) . This is 
an important reminder that the superfiuid non-axisymmetry 
and hence the gravitational wave strain depend on both the 
number of vortices that move and the degree of asymmetry 
in their positions within the star. We discuss this further in 
the context of real pulsar glitches in Sec. [S] 

The bottom panel of Fig. [T] is a spectrogram of |/ixi^, 
calculated by finding the power spectrum of \hx\^ in time 
windows of width At = 12.3. Since \hx\^ is calculated in 
the co-rotating frame, a peak at a; is absent from the ho- 
modyned power spectrum. Once again, bright features in 
the spectrogram (e.g. at i ~ 560) coincide with rotational 
glitches. The peak signal associated with each glitch occurs 
at the minimum preceding the glitch-induced spin up, for 
example, where the spikes in |/ix|'^ intersect the ijj(t) curve 
in the middle panel. In Fig. [21 we graph instantaneous power 
spectra from instants before, during, and after the glitches 
at i ~ 564 and i ~ 809 (black and grey curves respectively). 
All six curves exhibit a broadband burst signal, which in- 
creases with frequency, while the overall amplitude scales 
with |/ix 1^- The two solid curves exhibit peaks at frequency 
values 390 and 440, which are absent from the pre-glitch 
(dotted) and post-glitch (dashed) curves, suggesting that the 
peaks are associated with glitch-induced radiation. 
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The results presented in this section confirm that the 
motion of vortices through a grid of pinning sites results 
in a bursty gravitational wave signal. However, we empha- 
sise that a systematic study remains to be done of how the 
system size and range of pinning site strengths, amongst 
many other properties, change the gravitational wave sig- 
nal. Recent results show that once the ratio of pinning sites 
to vortices e xceeds unity , the statistics of spin-up events do 
not change (jWaJSzawski fc Melatos 2011' ) . However, when 
vortices outnumber pinning sites, pinning no longer dom- 
inates the spin-down dynamics of the condensate, result- 
ing in fewer and smaller glitches. An important unanswered 
physics question is, does a system containing ~ 10^^ vortices 
exhibit unpinning avalanches like systems with ^ 10^ vor - 
tices? Semi-analytic studies (jMelatos fc Warszawskil [2009|) 
suggest that avalanche dynamics are still possible in large 
systems, especially when unpinned vortic es can trigger fur- 
ther unpinnings in a domino-like effect (jWarszawski et al.l 
|2012| ). 



3 CURRENT QUADRUPOLE MOMENT 

Drawing inspir ation from the vortex d ynamics observed in 
recent studies (|Warszawski et al.ll2012f ) and the burst grav- 
itational wave signal computed in Fig. [T] we now develop 
an analytic formalism that facilitates calculation of a grav- 
itational wave signal from glitches involving realistic num- 
bers (> 10^) of vortice s. We follow the methodology intro- 
duced bv lThornj (jlOSOJ) and employed re cently in studies of 
f-mode and turbu l ence- d riven emission (IWasser manI 120081: 
Melatos fc Peraltal I2OI0I : ISennett et"all |a)10; Sid erv et all 
2OIOD . 

The far-field metric perturbation (wave strain) gener- 
ated by a superposition of current multipole moments S''" 
can be written as 
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in the transverse, traceless gauge, where t is the retarded 
time, D is the distance from source to observer, and Tj, ' ™ 
is the beam pattern, which is a function of the observer's 
orientation relative to the source. The physics of the source 
is housed in S'™. In general, hj^ also includes contribu- 
tions from mass multipoles j'™, which may dominate, but 
they are neglected in this paper, where we model the matter 
distribution inside a pulsar as incompressible and axisym- 
metric. 

The (i,m)-th mul tipole moment, S^^(t) (units: 
kgms"^), is given by (JThornel Il980l : iMelatos fc Peraltal 
l2010h 



S^ = ci I d xY'^r X ■ V X (pVs 
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where the integral is taken over the entire volume of the 
source, Vs and p are the fluid velocity and density respec- 
tively, V;^ is the complex conjugate of the scalar spherical 
harmonic 
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Pi" (cos 6) is the associated Legendre function, and one has 
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In our application, the leading term in the multipole ex- 
pansion is I = 2, for which only the m — +1 terms are 
nonzero [Piix) — — 3a;(l — x'^Y'^]. We make the simplify- 
ing assumption that the flow is purely azimuthaQ such that 
the mass current term x • V x (pVs) in Eq. Q reduces to 
z [V X (pvs)]^. It is convenient to switch to cylindrical coor- 
dinates {R' , (f)' , z' ) , such that, for a sphere of radius _Rs, we 
have 
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In Sec. |4] onwards, vortices are treated as discrete carriers of 
quantised circulation. It is then straight-forward to evaluate 
pVs and hence Eq. © as a sum over vortices. 

Denoting the inclination angle between the rotation axis 
of the star and the observer's line of sight by t, and the 
azimuth of the line of sight by (" (relative to some reference 
plane), the polarised components of the wave strain are 
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The angles t and ( enter through T^j, . For definiteness, we 
choose L — 7r/3 and ( = in most of the flgures and ap- 
plications to follow. In Sec. |4}[6l wave-strain estimates are 
presented mostly for /ix , as one has /i+ ~ /ix , except for cer- 
tain special (and unlikely) observer orientations. In App.lXl 
when calculating the stochastic background, we average over 
t G [0, tt] and set ( = without loss of generality. 



4 DISCRETE VORTEX MOTION 

A glitch occurs when many vortices simultaneously un- 
pin and move outward in response to a trigger. In this 
paper, wc remain agnostic about the natu re of the trig- 
ger, a notorious unsolved problem (lAnderso n & Itoh 1975|; 
Cheng et al.lll988l:|jahan-Mirill20()5l : lMelatos fc WarszawsM 
20091 : IWarszawski et al.ll2012l ). 



■^ iPeralta et al.l l l2006bl ) showed that this is inaccurate at high 
Reynolds numbers, when meridional circulation and Kolmogorov- 
like t urbulence set in (.Pcralta &i Mclatos.2009. : . Melatos fc Peraltal 
I2OIOI) . 
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Figure 3. Gravitational wave strain as a function of time in tlie cross polarisation hx (t) in units of Kq/t'^ for a single vortex with initial 
position Rf) = 0.1, which moves a distance Ar = 0.2, for angular velocity u)/{2-k) = 0.1, 0.5, 1.0 and 5.0 {top left, top right, bottom left 
and bottom right panel respectively). The vortex unpins at t = 1 and repins at t = 2 [bracketed by the vertical spikes in hx (*)]• 



Consider AA'v vortices that unpin from positions 
originally arranged in an evenly-spaced array [known as 
an 'Abrikosov lattice', which has hexagonal symmetry 
I Abrikosovl 119571 )]. A superfluid is irrotational; it only sus- 
tains vorticity at vortex singularities where one has 



V X V = ti5^ > [xx - xxv(t)] : 



(9) 



for a vortex whose axis runs parallel to the rotation axis and 
intersects the equatorial plane of the star at xxv(i). Here, 
K = h/m — 10~ m s~ is the quantum of circulation, where 
m is twice the neutron mass. The contribution to S^^ from a 
vortex at radius Rv{t) and azimuth (j}v{t) at time t in a fluid 
of uniform density p is therefore calculated by substituting 
Eq. (O into Eq. ^, with the result 
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Let us define dimensionless variables Rv = Rv{t)/Rs, 
t — t/r, (j)v = 4'v{i)/{2Ti), and Co = uir, where r denotes 
the glitch duration, and Rs is the stellar radius. The radial 
distance travelled by each vortex between unpinning and 
repinning is denoted by Ar = Ar/Rs (dimensionless). It is 
taken (arbitrarily) to be the same for all vortices and is a 
parameter of the model, to be evaluated in the future via 
Gross-Pitaevskii simulations. Hence we write Rv — Ro + 
d{t) and (pv = 4'o + Cot, where the initial values are chosen 
from (f>o € [0, 1) and Ro £ [0, 1 — Af] respectively. [Unless 
indicated, symbols used in plot labels are the dimensionless 
versions (tildes are omitted in plot labels).] In these units 



the single-vortex metric perturbation takes the form. 



, TT , 



ul^l^rn J. j-\ r7-,S2,21 -f^O -27ri<4o 



x^^ e 



'i?v(l-i?v)'' . (11) 



Non-glitch pulsar parameters are absorbed in the multiplica- 
tive constant 



G /512^y/' 4 



(12) 



except for Co. The glitch-related parameter r (vortex travel 
time) is kept separate. 

For a pinned vortex, with d(t) = or Af, the wave 
strain reduces to 



hjk {Ro,4'o,Cj,t) = -47r T^^ 



2rpB2,21 Kq -2Tri{4>o+u}i) 



xRo{l-R^y"Cj 



,2x3/2-2 



(13) 



Equation (I13|l is therefore the signal emitted by a vortex 
corotating with the stellar crust. Additionally, in the regime 
Cj ^ 1, where an unpinned vortex moves much further az- 
imuthally than radially during a glitch, Eq. (llSp is a good 
approximation to the wave strain resulting from the outward 
motion of a single unpinned vortex too. 
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Figure 4. Gravitational wave strain in tlie cross polarisation as a function of time hx{t) in units of Kq/t^ for a single vortex with 
uj/{2it) = 1.0 as a function of model parameters. Top left: Initial radius Rq = 0.1, 0.25, 0.5 and 0.75, with Ar = 0.2. Top right: Distance 
moved Af = 0.1, 0.2, 0.4 and 0.7, with Ro = 0.1. Bottom: Constant {solid curve) and parabolic {dotted curve) speed profile. All panels 
have irr = 1. 



5 SINGLE VORTEX SIGNAL 

In this section, we evaluate hx{t) for a single vortex as it 
unpins, moves, and repins. [The results for h+{t) are similar.] 
We investigate the effect on the wave form of changing ui, 
Rq, Ar, and d{t). The radial trajectory is described by 



d{i) 



It 





divv{i) 
Af 



t<U 



(14) 



t>U + l 



the vortex starts moving at t = ig. We experiment with both 
constant and parabolic speed profiles, «v(t), in what follows. 
The wave strain is calculated by substituting R^ — i?o + d(i) 
into Eq. ((TT1) . 

To begin with, we assume that the speed of the vortex 
whilst unpinned is constant, i.e. Vvit) = Af (tg ^ i ^ fg + 1 
and fg = 1). Figure [3] graphs the wave strain in the cross 
polarisation hx(i) in units of Kq/t^ [see Eq. (|12|) ] for 
u}/{2-k) = 0.1, 0.5, 1.0 and 5.0 (top left, top right, bottom left 
and bottom right respectively), corresponding to 0.1, 0.5, 1 
and 5 revolutions during the glitch. In all cases, tg = 1. 
The graphs extend one time unit beyond when the vor- 
tex repins. Even while pinned (t < 1 and f > 2), a single 
off-centre vortex produces an oscillatory gravitational wave 
signal with period 2ti/Cj, as its superfluid velocity field is 
not axisymmetric. Whilst the vortex remains pinned, the 
zero-to-peak amplitude of the wave strain, /imax, is pro- 
portional to up, as described by Eq. p3[l. While the vor- 
tex is moving, /imax changes at a rate 9(/imax/A"o)/9f = 



47r^cD2(l - Rlf^'^ \l - 2fv/(l - R^)] ■ We draw attention to 

the discontinuous gradient oi hx{t) at t = 1 and f = 2 (even 
when the vertical jumps are ignored), which occurs because 
the acceleration is artificially infinite instantaneously. Below 
we introduce a parabolic speed profile, which smoothes the 
acceleration and removes the spikes in hx (tg) and hx (tg + 1). 

The shape of the signal does not vary greatly with Rq- 
The top left panel of Fig. g] graphs hx for 0.1 ^ ^o sj 0.75, 
with fixed uj/{2n) = 1 and Af = 0.2. The amplitude does 
not vary monotonically with Ro, since the term pv, which 
appears in Eq. ([6]), is a function of the volume- weighted 
angular velocity (a single vortex does not generate rigid- 
body rotation). In fact, the amplitude peaks at Ro — 0.5 
followed by Ro — 0.25, 0.1 and 0.75 (in descending order). 
This ordering is preserved whilst the vortex is moving (1 < 
f < 2); however, for f > 2, the ordering becomes Ro = 
0.25, 0.1, 0.5, 0.75. Notably, the time at which the strain 
is a maximum for different Ro is not identical during the 
glitch. Once again, this feature arises because hx depends 
nonlinearly on R^; the sinusoidal signal resulting from the 
rotation of a pinned vortex is amplified or diminished as the 
vortex moves radially, resulting in a non-sinusoidal signal. 

The phase of the signal also varies with Af . We graph 
hx for different Af in the top right panel of Fig. [4] for 
u)/{2-k) — 1 and Ro = 0.1 fixed. For each value of Af , tra- 
versed in unit time, the in-motion wave strain scales pro- 
portional to Af . Once the vortex repins, the wave strain 



8 L. Warszawski and A. Melatos 



Quantity 


Value 


d)/{27r) 


1 


A(j> 


0.25 


Af 


0.1 


ACd/ui 


0.001 


AiVv 


500 



Table 1. Canonical glitch parameters adopted in the simulations 
in Sec. [B] 



amplitude depends on the new radial position _Ro + Af, as 
described by Eq. (|13|l . 

Finally, in the bottom panel of Fig. |4l we compare the 
wave strain for a vortex trajectory with a constant {solid 
curve; see above) and parabolic {dotted curve) speed profile, 



Mi) = 6Af (i ~ ig)[l - (i - ig)] 



(15) 



If we ignore the initial spike, a parabolic speed profile results 
in a higher wave strain, because the maximum v^ is 1.5 times 
higher than for constant speed. 

To summarise, the strongest gravitational wave signal 
from the motion of a single vortex is achieved when to and 
the maximum vortex speed are high, Af is large, and Rq — 
0.5. 



6 VORTEX AVALANCHE SIGNAL 

6.1 Unpinning geometry 

Current glitch detection methods based on radio telescope 
data are sensitive to glitches involving the simultaneous un- 
pinning and outward motion of at least AA^v ~ 10^^ vortices. 
The following approximate formula relates AA'^v to the frac- 
tional frequency jump Auj/u observed: 



Aoj 



Af 



AJVy Js 

iVv 7c 



(16) 



It involves the ratio of the superfiuid and crust mo- 
ments of inertia, Is/Ic (which depends on what fraction 
of vortices are pinned; Ib/ Ic ~ 10~^) (JAlpar et al.l 1 1989 : 
lEpstein fc Van Ripen Il992l '). and the total number of vor- 
tices A^v « 2iiwRI/k « lO^'^ [aj/(l rad s"^)]0 Hence 
AA'v = 10^^ corresponds to Auj/lj ~ 10~^^, typically the 
small es t glitch size o b served to date (Ijanssen fc Stapperd 
20061: iMelatos et al.l l2008l : IChukwude fc Uramal l2oTa 



Espinoza et al.ll201ll 'l 



The geometry and chronology of a vortex avalanche 
are vital, yet regrettably unknown, ingredients of the grav- 
itational wave calculation. Do all vortices move simul- 
taneously, from evenly distributed positions throughout 
the star? Evenly distributed pinning sites with a broad 
range of pinning potentials represent one possible con- 
figuration that would give rise to this scenario l|jonesl 



■^ This approximation assumes that each vortex is responsible for 
an equal fraction of the total angular momentum of the superfiuid, 
which is not accurate when the vortices are spread evenly from 
the centre to the edge of the star. 



Il99ll : iMelatos fc WarszawskH |2009| '1. Or, do vortices un- 
pin in a cascade restricted to a wedge with in the star 
(|Alpar et al.lll98ll :IWarszawski fc Melatos|[200i ')? A branch- 
ing process, where an unpinned vortex also preferentially 
unpins its nearest neighbours, leading to a branching tree 
of activity, would produce th i s out come. We adopt the 
terms 'creep-like' (lAlpar et al.l 19891 ) and 'avalanche- like' 
(see IWarszawski fc Melatoa [2OO3) to describe the for- 
mer and latter possibilities respectively. During creep-like 
glitches, the non-axisymmetric change in velocity field due 
to one vortex largely cancels a similar change from a vortex 
on the opposite side of the star. Hence we have (^x) =0 (en- 
semble average) and the typical strength of the signal in any 
particular realisation is of order the dispersion (/ix)^ • Dur- 
ing avalanche-like glitches, cancellation is incomplete and 
(hx) is not zero in general. 

In what follows, we calculate only the contribution to 
the wave strain from vortices that unpin and move during 
a glitch. In reality, even the static, pinned Abrikosov vortex 
lattice is not perfectly smooth and axisymmetric: vorticity 
is carried by point-like vortices rather than in a smooth 
field. We ignore the small contribution to the wave strain 
from pinned vortices. In addition, we neglect the contribu- 
tion from persistent asymmetries in the vortex distribution 
near the stellar c ore due to incomplet e, stratification-limited 
Ekman spin up (jBennett et al.11201 0'). 

We perform Monte-Carlo calculations by choosing the 
initial positions of unpinned vortices at random to inves- 
tigate the dispersion in hx as well as its mean. The left 
panel of Fig. [5] shows the initial vortex configuration for 
A'v = 10*. Vortices (represented as dots) are initially placed 
in an hexagonal Abrikosov lattice, inside a circle of radius 
1 — Af . The centre panel of Fig. [5] shows the vortex po- 
sitions after a creep-like glitch. The right panel shows the 
positions following an avalanche-like glitch of the same size. 
The initial positions (_Ro and <^o) of vortices that unpin are 
drawn from the following probability distribution functions 
(PDFs): 



p{Ro) 



2Ro 
1 



(17) 



q{M = -^ [H{<t)o + A<^o/2) - H{4,o - A0o/2)] , (18) 
A0O 

where H{) is the Heaviside step function, and we choose 
00 = as the bisector of the avalanche without loss of gen- 
erality. For creep-like and avalanche-like glitches, we have 
A^o = 1 and A4>q < 1 respectively. We adopt the parabolic 
speed profile from Eq. (|15p . so as to avoid discontinuities 
leading to artificial spikes in the wave strain. A parabolic 
profile is also a close approximation to the spiral motion of 
a vortex as it repins, observed in the quantum mechanical 
simulations described in Sec. |51 



6.2 Wave strain 

We now calculate the gravitational wave signal from a vortex 
avalanche. We evaluate hx{t) in two ways: (1) Monte-Carlo 
simulations, in which AAV, Ro and (jio are drawn from the 
avalanche size PDF, g{Au}/ui), p{Rq) and q{cf>o) respectively, 
using Eq. (|16p with Is/Ic = 1; and (2) marginalising over 
the PDFs in Eq. (fT7|l and ([TSJI to obtain the moments of 
hx{t) for fixed Auj/uj. Since the contribution to hx from 
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Figure 5. Left: Initial vortex positions for 10* vortices in an equatorial cross-section of a neutron star. Vortices are confined to radii 
_R < 1 — Af so that unpinned vortices do not leave the star. Centre: Final vortex positions after a creep-like glitch. Right: Final vortex 
positions after an avalanchc-likc glitch with opening angle A(f> = 1/8. 
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Figure 6. Monte-Carlo results for the wave strain in the cross polarisation as a function of time hx (i) for vortex avalanches with opening 
angles Ac^o = 1-0, 0.95 and 0.35 {top left, top right and bottom panels respectively). Solid grey curves represent hxit) for 10 different 
glitch realisations. Dashed grey curves display the mean plus/minus one standard deviation of 20 (only 10 are plotted) Monte-Carlo 
realisations. Dashed black curves are theoretical predictions of the mean plus/minus one standard deviation given in Sec. 16.21 Simulation 
parameters: TVv = 5 X 10"*, i:j/(27r) = 1.0, Af = 0.1. 



each vortex is drawn from the same underlying PDF, and 
these contributions are summed together, the statistics of 
the aggregate /ix it) obey the central limit theorem (AN^ 2> 
1). That is, at each time t, hx (t) is a Gaussian variable with 
mean AN^fJ-i and variance ANvO-f, where /ii and af are the 
mean and variance respectively of the contribution from a 
single vortex, given by 



a?(i) 



2Ko 
r^A(f)Q Jo 
(hi) - fil 



dRo Ro 



A0/2 



A</i/2 



d(j)oh 



jk 



(19) 
(20) 



(hi) it) = 



K0V2 
T^A4,Q 



dRo Ro 



A0/2 _ 2 

dcpolkj^m) 

-A0/2 ' ' 



Here, hjl[ is calculated from Eq. pip . For A(f>o = 1, 
Eq. (|19p implies fii = 0. However, any particular realisa- 
tion of the initial vortex positions is not symmetrically dis- 
tributed around the star, implying af > 0. 

In Fig. [6] we graph hx{i) (in units of Ko/t^) from 10 
glitch realisations {grey curves), for Acj) = 1, 0.95 and 0.35 
[top left, top right and bottom panels respectively), overplot- 
ted with curves representing one standard deviation above 
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Figure 7. Gravitational wave strain statistics from Monte-Carlo simulations in units of KqJt^. Top left: Maximum {^crosses) and 
time-averaged (asterisks) (hx) (in units of Kq/t'^) as a function of avalanche opening angle Ac^q- The solid and dotted curves are 
the corresponding analytic predictions. Top right: Maximum (crosses) and time-averaged (asterisks) standard deviation [crjj = ((/i^) ~ 
C'x)^)^ ] as a function of avalanche opening angle A0o- The solid and dotted curves are the corresponding analytic predictions. Bottom: 
Maximum (crosses) and time-averaged (asterisks) standard deviation as a function of number of vortices involved in the avalanche, AA'^v . 
The solid and dotted curves are the corresponding analytic predictions. 



and below the mean (dashed grey curves). Unless other- 
wise stated, we adopt the set of canonical glitch parame- 
ters in Table [T] We compare these curves to the analytic 
prediction AN^fi-i{t) ± y/AN^ai(t), given by Eq. (fT9)) and 
H20[) (dashed black curves). Both methods of calculation as- 
sume that all AA'v vortices move simultaneously; in a real 
avalanche, vortices unpin in a domino chain, but the domino 
time-scale is rapid compared to the time separation between 
glitches. Even for 20 realisations, the standard deviation of 
the Monte-Carlo results, (tmc, deviates from aiVAA'v by 
less than 1% for A0 — 0.35. For A(/!> = 1, the realisations 
yield a mean signal that oscillates about zero [Eq. (|19|l gives 
precisely zero; our simulations involve a finite number of 
vortices, resulting in residual dispersion]. 

The position of the strain maximum, or phase of the 
signal, varies by up to half a revolution from glitch to 
glitch. The phase shift appears for the same reason dis- 
cussed in Sec. [S] the signal amplitude from each vortex in 
the avalanche is a non-monotonic function of its radial po- 
sition. For a given AA'v, avalanches that are restricted to 
a thin wedge of the star are more likely than avalanches in 
larger regions to involve vortices from all radii, effectively 
smoothing over this effect; for A(j} <^ 1, the phase difference 
between different realisations is negligible. 

In order to systematically parametrise the gravitational 
wave signal, we define four quantities ((■) denotes an ensem- 



ble average over realisations and ah is the standard deviation 
of the glitch signal from many vortices). 

(1) The time-averaged mean, (hx) = t~^ J^ dt{hx){t). 

(2) The maximum amplitude of the mean {/ix)max ~ 
max{{hx)(t)|te [0,r]}. 

(3) The time-averaged standard deviation, (ah) ~ 
r-' j; ah(t)dt. 

(4) The maximum amplitude of the standard deviation, 
o-max = max{ah|t G [0,r]}. 

The top panels of Fig. [7] graph these parameters as a 
function of A0. For (hx) and {/ix)max, Eq. (|19p conforms ex- 
actly to the Monte-Carlo output (see left panel). The signal 
strength is a decreasing function of the avalanche opening 
angle: as Aff) increases, it is more likely that the signal from 
the motion of one vortex interferes destructively with that 
of another vortex on the other side of the star. The top right 
panel shows that, for A0 > 0.3, (hx) decreases linearly with 
A0. Both (ah) and a^iax are increasing functions of A0 for 
A(jj < 0.5. For Ac/) > 0.5, (ah) and amax saturate, since the 
avalanche spans more than half the star, and hence the mo- 
tion of vortices on one side of the star start to 'cancel' with 
motion of vortices on the opposite side. 

The bottom panel of Fig. [7] graphs (ah) and cTmax as 
functions of AA'^v. As predicted by the central limit theorem, 
both quantities increase with the square root of AA^v. For 
glitches involving a significant fraction of the vortices, this 
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top panels of Fig. [S] we graph the mean {left panel) and stan- 
dard deviation {right panel) of the wave strain as a function 
of glitch size, expressed in terms of AA'v, for avalanche-like 
glitches with A^ — 0.25. Both measures of the mean {top 
left panel) introduced in Sec. l6.2l are linearly increasing func- 
tions of ANv The Monte-Carlo {crosses and asterisks) and 
analytic {solid and dotted curves) calculations agree with re- 
spect to the means. On the other hand, analytic calculations 
of (ah) are up to double the Monte-Carlo values, and cTmax 
does not agree well with Eq. (|20p . although all calculations 
respect the scaling (Tmax oc y/AN^. The same scaling is ob- 
served for creep-like glitches, as shown in the bottom panel 
of Fig. [3 
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Figure 9. Monte-Carlo results for the wave strain in the cross 
polarisation as a function of time hx{t) for angular velocity 
ui/{27t) = 0.1, 1.0 and 5.0 {top, middle and bottom respectively). 
Solid grey curves represent hy,{t) for different glitch realisations. 
Dotted grey and dashed grey curves are the mean and the mean 
plus/minus one standard deviation of 20 (only 10 are plotted) 
Monte-Carlo realisations. Dotted black and dashed black curves 
are the corresponding analytical predictions. Simulation parame- 
ters: Afv = 5 X 10*, A^o = 0.25, Af = 0.1. 



trend breaks down, since most of the vortices unpin during 
the avalanche, not just a small random sample. 

In the following sections, we consider avalanches with 
A(f> = 0.25, and explore the change in the mean and stan- 
dard deviation of the signal as a function of pulsar angular 
velocity, glitch size, and vortex travel distance. 



6.3 Glitch size 

Both avalanche and creep-like glitches generate wave strains 
that scale linearly with glitch size (and hence AA'^v). In the 



6.4 Vortex travel distance 

Glitch size is the product of the number of vortices that un- 
pin and the distance they traverse, according to Eq. (|16p . 
In other words, for a given glitch size, either many vortices 
move a short distance, or few vortices move a large distance. 
In this section we investigate how the strength of the gravi- 
tational wave signal changes with the distance travelled by 
each vortex, Af, to cause a glitch of fixed size. 

The centre panels of Fig. |S] demonstrate that the mean 
gravitational wave strength is a decreasing function of Af. 
This result can be understood by noting that increasing the 
radial distance traversed by each vortex does not make the 
system more or less asymmetric. Therefore, as Af increases, 
and AN^ decreases in inverse proportion to keep Auj /uj con- 
stant, we expect the mean and standard deviation to also 
decrease. This is confirmed by the Monte-Carlo and analytic 
results in Fig. [S] In calculations not graphed here, we also 
find that this correlation persists for creep-like glitches. 



6.5 Pulsar angular velocity 

The mean and dispersion of the wave strain increase with 
oj; (hx), (/ix)max, {(^h) , aud amax are plotted versus Cj in 
the bottom panels of Fig. [8] The quantities are quadratic 
functions of u, as for a single vortex, e.g. Eq. (|13|l in Sec. U) 

Figure |5] demonstrates that, for glitches in which az- 
imuthal vortex motion is faster than radial motion, the grav- 
itational wave signal oscillates with the same angular veloc- 
ity as the pulsar and diminishes with time. For Cj/{2tt) < 1, 
the signal executes only a fraction of a full oscillation. The 
variation in the signal results from both azimuthal vortex 
motion and the changing parabolic vortex speed profile. For 
example, in the top panel of Fig. [O] the signal executes ap- 
proximately half an oscillation, and the vortex circulates 
one tenth of the way around the star [uj/{2n) — 0.1]. The 
dependence on uj is most obvious in the bottom panel of 
Fig. |9l which graphs hx{t) for Cj/{2tt) = 5.0 [solid grey 
curves are Monte-Carlo realisations, dashed curves are the 
mean plus/minus one standard deviation for the Monte- 
Carlo {grey) and analytic {black) calculations]. Over the 
lifetime of the glitch, the amplitude of {hx){t) decreases by 
20%. 

From Sec. [5] we know that the signal from a single vortex 
is a function of its radial position. As a glitch progresses, the 
mean (taken over all vortices) radial position increases from 
its initial value (Rv) = 2/3. Figure |4] demonstrates that. 
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Figure 8. Gravitational wave strain statistics from Monte-Carlo simulations of many-vortex glitches in an individual pulsar. Maximum 
{crosses) and time-averaged (asterisks) {hx) {left) and ctj^ {right). The solid and dotted curves are the corresponding analytic results. 
The mean and standard deviation are given in units of Kq /r ^ . Top : Variation with the number of vortices that move during the glitch, 
AA''v. Centre: Variation with distance travelled by each vortex, Ar. Bottom: Variation with angular velocity, ui. 



for 7?v > 0.5, the maximum signal strength is a decreas- 
ing function of i?v, and hence we expect the signal from 
an avalanche to weaken with time. It is important to note 
that for a real glitch this effect is expected to be negligi- 
ble (Af ^ (Ro))- As in previous sections, the number of 
oscillations of {hx}{t) during the avalanche is not precisely 
uj /{2n), owing to changes in (Rv) with time. 



vortices still travel across only across a very small fraction of 
the star), the azimuthal motion of the vortices is negligible 
compared to their radial motion; this scenario arises when 
the vortex travel time is much shorter than the neutron star 
rotation period. In this case, the zero-to-peak amplitude of 
hx is obtained by setting dR^/dt ~ Af/r and expanding 
Eq. (|11|) to first order in Af and is given by 



6.6 Amplitude scaling 

As a guide to how the signal from a single glitch scales with 
the physical parameters associated with the glitch, we pro- 
vide a back-of-the-envelope estimate for the maximum wave 
strain. There are two regimes, determined by the relative 
size of UJ and Ar. For cD <^ Af, but keeping Af ^ 1 (i.e. the 
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3^Af{RofAN. 
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Ikpc/ V 

Alo/lj 



10 ms 



10- 



VlO^rads" 



(22) 



(23) 



Gravitational-waves bursts from vortex avalanches 13 



D and lj are the distance to the source and its angular ve- 
locity respectively, and they are pulsar, but not glitch, de- 
pendent. T and Ar are the vortex travel time and distance 
respectively. We note that both r and Ar may vary be- 
tween glitches, and even between individual vortices moving 
within a single glitch. However, in general, the faster and fur- 
ther vortices move, the stronger the gravitational wave sig- 
nal. Further studies using first-p rinciples simulations, suc h 
as those presented in Sec. \7\ and IWarszawski et al.l ([2013), 
and an understanding of how vortex behaviour scales with 
system size, are needed to fix these parameters accurately. 
Alternatively, for uj ^ Ar, but to still small, the az- 
imuthal vortex motion is a larger contributor to the wave 
strain than is the radial motion. The signal amplitude is ob- 
tained by expanding Eq. (|11[) to first order in u, and scales 



{f^X ,inax/ 



^{27Vi:jf{Ro}AN. 



= 10 
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Ar 



10- 
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/Aw/ojX / IjJ 

V 10-^ ) VlO^rads-^ 



(24) 



(25) 



In this regime, corotation of pinned vortices with the crust 
produces a non-zero signal, provided the superfluid flow is 
non-axisymmetric, whose strength is a f unction of the non - 
uniformity of the vortex configuration (| WassermanI l200a ) . 
e.g. pinning at grain boundaries (inhomogeneous) versus dis- 
locations and lattice sites (homogeneous). Here, it is impor- 
tant to note that Eq. (|25p estimates the strength of the 
gravitational wave strain from azimuthal vortex motion on 
the rotation time scale, even though such motion is ongoing, 
even when the neutron star is not glitching. To avoid the con- 
clusion that there is a strong, persistent signal from rotation 
of pinned vortices, we must assume that non-axisymmetry 
in the vortex lattice is much more pronounced when a large 
number of vortices are not pinned (i.e. during a glitch) than 
during inter- glitch intervals. 



6.7 Example: the 2006 glitch in Vela 

We can use Eq. (|23p and (|25p to constrain the glitch dura- 
tion for the 2006 Vela gl itch. Non-detection of gravitational 
waves from this glitch (JAbadie et al.l 1201 ll ) put an upper 
limit of 10^"^" on the wave strain. For uj <C Af, with Vela- 
like parameters {D « 300 pc, uj « 71rads~^, Auj/uj ~ 10~^ 
and Ar = 10~^m), the lower bound on the glitch dura- 
tion is then r > 10"* ms. Of course, the independence of r, 
Ar and Auj/uj cannot be taken for granted. For example, 
if T is the time taken for a vortex to travel a distance Ar, 
which in turn depends on the superfluid-crust velocity lag 
(through the vortex velocity), then we cannot hold Ar, and 
even Alj/uu, constant whilst varying r. We also remind the 
reader that underpinning these estimates is the conservative 
assumption that the pinned vortex distribution is symmet- 
ric, apart from those vortices participating in the avalanche. 
Otherwise there would be a persistent gravitational wave 
signal described by Eq. (P5|) . 

We also note that the uj ^ Af approximation returns 
a wave strain that is much stronger than the persistent sig- 
nal (h ~ 10~^^) arising from turbulence in a differentially 



rotating Crab- like neutron star IjMelatos fc Peraltall2010l l. 
The predicted wave strain resulting from the same physical 
mechanism operating in a nearby {D ~ 0.01 kpc), rapidly 
rotating (il ~ 3 x lO^rads"^) neutron star is comparable 
to our UJ ^ Af scenario [Eq. ()25|l ] for a glitching pulsar 
(|Melatos fc Peralta|[201ol '). 



7 DISCUSSION 

This paper estimates the strength of the gravitational wave 
burst arising from individual pulsar glitches. Quantum me- 
chanical simulations of rotating, pinned superfiuids demon- 
strate that a broadband burst of gravitational radiation is 
emitted in the current quadrupole channel by the spasmodic, 
non-axisymmetric rearrangement of superfluid vortices dur- 
ing a glitch. Generalising the simulation results, which are 
restricted to systems with ~ 10^ vortices, we employ a dis- 
crete model of vortex motion to calculate analytically the 
wave strain from vortex avalanches involving realistic num- 
bers of vortices (~ 10^^). 

When a single vortex unpins and moves, the wave strain 
peaks when the angular velocity and the distance travelled 
are large, the initial position is half the stellar radius, and 
the vortex speed is high. Monte-Carlo simulations agree well 
with analytic calculations. The wave strain depends strongly 
on the unpinning geometry. If vortices unpin at random 
throughout the star (vortex creep), the signal is smaller 
than when unpinning is confined to a thin wedge (vortex 
avalanche); the maximum wave strain from a glitch of fixed 
size varies inversely with Acp. In fact, for A<^ = 1, the mean 
wave strain is formally zero. 

Quantitatively, the maximum strain is inversely propor- 
tional to Af , quadratic in uj, and linearly proportional to the 
number of vortices that unpin and move. The standard devi- 
ation of the maximum strain depends on these physical pa- 
rameters in the same way. We note, however, that for a glitch 
of a given size, Af and AA'v are not independent variables. 
That is, for fixed Af , the number of vortices that unpin and 
move during a glitch of size s is AA^v ~ sAv/c/(/sAf). We 
also find that the increase in average radial vortex position 
as vortices move outward during a glitch damps the ampli- 
tude of the glitch gravitational wave signal. This is most ap- 
parent for glitches with large uu; the glitch signal oscillates, 
and its zero-to-peak amplitude decreases with time. How- 
ever, in real glitches, the change in average vortex position 
during a glitch is so small that this effect is negligible. 

The wave strain for a glitch in which a vortex travels 
further azimuthally than radially (ui S> 1) has a zero-to- 
peak amplitude h ~ 10~^*. For a glitch dominated by radial 
motion (uj -^ 1), the zero-to-peak amplitude is /i ~ 10"^*. 
These estimates are based on a glitch of size 10~^. Conve- 
nient scalings of wave strain with glitch size, stellar angu- 
lar velocity, vortex travel distance, and glitch duration are 
presented in Eq. (|23|) and (|25|) respectively. We remind the 
reader that we conservatively assume that vortices that do 
not move during a glitch do not contribute to the gravita- 
tional wave signal. Therefore, when calculating the signal 
from oi-dominated vortex motion, we implicitly assume that 
vortices involved in a glitch are asymmetrically distributed. 

Unlike the individual glitch burst signal from a nearby 
pulsar, which is eminently detectable, the contribution to 
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the stochastic gravitational wave background is arguably 
weak. Nevertheless, for completeness, in App. |XJ we ex- 
plain how the theoretical machinery in Sec. |3]-|6] can be 
applied to calculate the contribution to the background 
rigorously. Estimates of the strength and spectrum of a 
stochastic background of gravitational waves resulting from 
sources of cosmological and astrophysical origin are essen- 
tial to determining the det ectability of bu rst events and con- 
tinuous wave sources f] Allen fc Romandll999 : Ferrari et al.l 
1 19991 : iRegimbau fc de Freitas Pachecdl2006bl i. We find that 
power-law-distributed glitch sizes (power-law index —3/2) 
from a Galactic population of neutron stars contribute to the 
stochastic background with a signal-to-noise ratio ~ 10~^ 
(Einstein Telescope). Importantly, however, the inclusion of 
extra-Galactic pulsars must boost the signal substantially, 
something that merits careful consideration in a future pa- 
per. 

The calculations presented here are premised on sev- 
eral assumptions about the physics of glitches and the glitch 
behaviour of the local neutron star population. Most im- 
portantly, our calculation assumes that glitches result from 
superfluid vortex unpinning avalan ches. Alternative glitch 
theories, such a s fluid instabilities (lAndersson et al.l |2004| ) 
and starquakes (JLink fc Epsteinlll996l : lNegill201(J ). generate 
different signal strengths and spectra. Previous estimates of 
the maximum achievable signal-to-noise ratio for glitch de- 
tections from the energy deposited in /-modes give signal- 
to-no ise ratios 0.4 < S/N < 3 x 10^ JAndersson fc Comeil 
I2OOII '). The range in S/N is given for a EURO detector 
including (lower bound) and ignoring (upper bound) shot 
noise. Larger S/N (up to ~ 7.4 x 10^) are predicted for in- 
stability modes in which the pro ton and neutron fluid s in the 
stellar core are counter-moving. lAbadie et al.l (|201in placed 
Bayesian 90% confidence upper limits on the peak ampli- 
tude from the 2006 Vela glitch of 1.4 x lO"^'', corresponding 
to a maximum energy released through oscillations of the 
fundamental quadrupole mode of 1.3 x 10'^^ erg. 

Even within the vortex unpinning paradigm, changes to 
the unpinning physics may alter the estimated gravitational 
wave signal strength significantly. If, as suggested in the pul- 
sar g litch literature (jMelatos et al.l 120081 : 1 Warszawski et al.l 
120121 ). glitches result from a domino effect, then contribu- 
tions from individual vortices should be summed consecu- 
tively and coherently, rather than simultaneously. This is 
likely to lead to a greater wave strain, as the motion of one 
vortex is less likely to cancel geometrically with another. 

In conclusion, the gravitational wave signal emitted dur- 
ing the spin-up phase of individual glitches resulting from 
superfluid vortex avalanches is not detectable by current de- 
tectors, but it is close as Eq. (|25p indicates. Once a detec- 
tion is made, other informative experiments are possible, 
e.g. searching for all-sky correlations between the -|- and x 
polarisations. These matters are deferred to future work. 



APPENDIX A: STOCHASTIC BACKGROUND 

Gravitational waves from neutron star glitches contribute 
to a stochastic gravitational wave back ground composed o f 
many unresolved astrophysical sources ( Abbott et al.ll2009f ) . 



far (jManchester et al.l 120051 '). Making some simple assump- 
tions about the glitch rate of the neutron star population 
as a whole (non-pulsating objects should emit gravitational 
radiation when they glitch too), we estimate the strength 
of the glitch background. In this first attempt, we ignore 
the long-term spin-down evolution of pulsar spins during 
the lifetime of the Galaxy. We also restrict attention to 
sources in the Milky Way, as we know nothing about the 
glitch properties of extra-Galactic neutron stars. With the 
latter restriction, we find that the Milky Way glitch back- 
ground is weaker than other LIGO backgrounds. However, 
the theoretical machinery to calculate the signal, which flows 
directly from the analysis in Sec. [2][6l is described here for 
completeness. Simple order-of-magnitude estimates suggest 
that the background may be significant when extra-Galactic 
contributions are included, a topic for a future paper. 



Al Total glitch rate 

We estimate the total rate of glitches, O, by assuming that 
glitches are manifestations of a self-organised branchi ng pro- 
cess, char acterised by a power-law size distribution (jHarrig 
[l989; Dal v fc Porporatd 2007) 



g{s) = (-l/2)(6 
where Smin ^ 



-1/2 



-1/2n 



^1^-3/2 



(Al) 



^ s = Aijj/ctj ^ s+ is the range of frac- 
tional glitch sizes. This distribution is supported quali- 
tative W_b^;_observational data for several individual pul- 
sars (jMelatos et al.l 120081 ) . although it can be shown by a 
Kolmogorov-Smirnov test on aggregate data that g{s) oc 
s^'^'^ is not universal (Melatos et al. 2008). Across the pop- 
ulation, one observes s+ ~ 10"* and Smin ~ 10~^^ (the ob- 
served minimum glitch size) respectively. The physical lower 
bound on glitch size, s_, is given by the fractional change in 
angular velocity induced by the outward radial motion of a 
single vortex by a distance characterised by the mean inter- 
vortex separation Ar ~ [k/{2uj)]^'^ . Since the total number 
of vortices in a pulsar is A^v = 2itujR'^/k, we obtain 



(JL) 

\2ujJ 



1/2 



2nuiR^ Ic 
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(A2) 
(A3) 



The Milky Way is populated by ~ 10 neutron stars, of 
which ~ 2 X 10"^ have been observed electromagnetically so 



To proceed, we ask the question, what fraction of 
glitches are within the observable range? The integrated 
probability between Smin and s+ is 

r dsg{s) = (s-[f -.;i/^)/(s;^/^-s:^/^) ~ lo-^ .(A4) 

Hence only one in 10^ glitches is observed by radio tele- 
scopes at current sensitivities and for typical observational 
duty cycles. We now make the simplifying assumption that 
the glitch rate in each individual pulsar. A, falls off ex- 
ponentially with characteristic age [A(rc) oc ex p(— Tc/Tp), 
where Tp is some referen ce characteristic age] (JLvne et al.l 
l200d : iMelatos et al.lliooi '). Tc = uj/{2dj), and that neutron 
stars are being born at a constant rate, leading to a fiat 
distribution of neutron star ages in the galaxy, p(T'c) = 
{Tmax — rmin)~^- A coroUary of this first assumption is that 
the Grab pulsar (Tc ~ 10^ yr) glitches approximately three 
times more often than the Vela pulsar (Tc « 10^ yr), which 
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agrees roughly with observations IjMcKenna fc Lvnelll99CI : 
iMelatos et aill2008l 'l. The PDF of gUtch rates is then 



P(A) = -pr (Tnmx — Tmin 



(A5) 



Normahsation gives A « 10'' yr~^, where Tmin ~ W^ yr and 
Tmax ~ 10^" yr are the minimum and maximum age of local 
neutron stars respectively. The mean glitching rate is then 

(A) = / dX Ap(A) (A6) 



J- min 
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-^ min 
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which translates into (A) ~ 10~^'^s~^ for an individual pul- 
sar. In other words, the mean is dominated by old objects. 
Given a population of TVp = 10® pulsars in the Milky Way, 
the total glitching rate is then 



e = loM^ 



(A) 



10- 



(A8) 



[the factor of 10^ comes from Eq. (|A4[) . which accounts for 
unobserved small glitches]. In the next section we see that 
the relative size of 6 and r determines the smoothness of 
the gravitational wave background. 

For simplicity, and due to a lack of information, we ex- 
clude glitches in extragalactic pulsars from the above dis- 
cussion. However, we show below that the energy density in 
gravitational waves is proportional to the aggregate glitch 
rate in the cosmic volume under consideration, Q, divided by 
(D)^, the square of the average distance between source and 
detector. One typically has Q oc (D)^ for isotropically dis- 
tributed sources, suggesting that the contribution from ex- 
tragalactic pulsars actually dominates. We defer this impor- 
tant calculation to future work. To do it properly, one must 
account for the redshift-dependent source formation rate 
and in put the statistical properties of extra-ga l actic neutron 
star s ( Rcgim bau fc de Freitas Pachecd l2006al : I Howell et al.l 
[20rit : .Marassi et al.ll2o"lol '). 



A2 Temporal structure 

The character of the stochastic gravitational wave back- 
ground from pulsar glitches depends on the average number 
of concurrent glitches, given by Or. In the regime 0r <^ 1, 
glitches produce a 'shot noise' background of distinguishable 
burst events. In Sec. l6.6l we place upper and lower bounds on 
the strength of such distinguishable events. In the opposite 
regime, 0r ^ 1, glitches generate a continuous (i.e. never 
absent) background, whose mean and variance we estimate 
in App. IA3I The intermed iate regime &t ~ 1 is often de- 
scribe d as 'popcorn' noise IjRegimbau fc de Freitas Pachecd 
l2006al ). 

Let us construct the functional form of the background 
in a time interval of length T, during which ST bursts arrive 
at Earth. The wave strain during this interval is 

hjifit) = hi(t - ti) + hiit ~ fa) + ... + herit - ter) , (A9) 

where t„ is the start time of the n-th glitch, and hn(t — t„) 
is the gravitational wave signal from the n-th glitch. We 
make the simplifying assumption that the glitch duration, 
r, is the same for all glitches, and that the glitch profiles 



hi{x) = ... — hQT{x) have the same shape but different 
amplitudes. Continuous and shot noise arise from situations 
with {f„— f„_i) <^ T and {tn — tn-i) ^ r respectively, where 
the averages are taken over many inter-glitch intervals. 

For each glitch, we select values for s and D from proba- 
bility density functions describing the distributions of glitch 
size and distance from Earth, g{s) and A(D). The latter 
PDF is 



HD) = Q 



D + Di 



exp 



D-Dq 

Dq+Dj 



(AlO) 



where Q is a normalisation constant, a = 1.64 ± 0.11, 
b = 4.01 ± 0.24, Di = 0.55 ± 0.10 kpc, Dq is the Earth- 
Sun dis tance and D is measured from the centre of the 
galaxy (|Faucher-Giguere fc Kaspil |200a). The inter-glitch 
time is chosen from an exponential PDF with mean rate O 
JMelatos et al.ll2008l ). We use Eq. (O to evaluate /in(i-f„). 

Figure [M\ graphs hx{i) for Gr = 1 and 10^ {left and 
right respectively) and Lb/{2-n) — 1.0. For 0r = 1, each 
glitch is distinguishable. On average, ©r glitches are stacked 
at any given time. Glitch sizes are drawn from a power law 
that covers 22 decades. Feasible Monte Carlo simulations do 
not sample this distribution comprehensively, so that mean 
sampled glitch size is many orders of magnitude smaller than 
Jdssg{s). However, Monte-Carlo simulations used here to 
check the analytic theory, not to deliver actual numbers, 
which is deferred to a future paper. For testing purposes, a 
compressed s range is ample. 

In the left panel of Fig. lA2l we plot the mean wave strain 
as a function of Or for a superposition of glitches whose sizes 
are in the range 5 x 10~^ to 1 x 10~^. Again, the reduced 
range ensures that the Monte-Carlo simulation comprehen- 
sively samples g{s). A least-squares power-law fit returns a 
power-law index of 0.94. The mean standard deviation {(7^), 
graphed in the right panel of Fig. IA2I is also well described 
by a power law with index 0.64, roughly consistent with 
{ah) oc VAA^v 



A3 Gravitational wave energy density 

We now use the analytic expressions derived in Sec. |6] to 
calculate the glitch contribution to the dimensionless, cos- 
mological, energy density parameter figw. We choose this 
quantity to make contact with other studies of stochas- 
tic gravitational-wave backgrounds, includi ng the recent 
LIGO bound on theories of cosmic strings (jAbbott et al.l 
|200 9|). The calcula tion is relevant to the regime Or ^ 1 
(; Ho well et al.l 1201 1 | ). in which the bac kgro und is con t inuou s 
in time. Following Ferrar i et all ([l99i) and lMaggior3 (|200Ci ). 
we write Slgw as a function of wave frequency, v, by relating 
it to the spectral energy density, d'^Eg„/{dudS), through the 
expression 



a 



gw(j^) 



I/O d^E^ 



C^pc 



dvdS 



(All) 



where per = SHq/{8tvG) is the critical density to close 
the universe, Hq is Hubble's constant, is the total event 
rate from Sec. lAll dS is the area element on the sky, and 
-Egw {i^, n) is the total energy radiated in the direction n at 
frequency u. 

We begin the task of calculating d?E^^/{dvdS) by writ - 



we Degm tne task ot calculating a £jf,„/(auaoj Dy writ - 
ing the energy flux from a source as IjKokkotas et al.ll200ll ) 
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Figure Al. Top: Gravitational wave strain in tlie cross polarisation hx{i) in units of Kq/t^ as a function of time from glitches from 
multiple pulsars. The time interval is sufficient to bracket five consecutive glitches. Glitch sizes are drawn from a power-law distribution 
with minimum and maximum s^i^ = 5 X 10~^ and Smax = 1 X 10~® respectively. The time between glitches is drawn from an exponential 
distribution with mean rate = 1 (left) and = 10^ (right). Parameters: Af = 10~^, t = 0.01, ui/(2tt) = 1. 
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Figure A2. Mean (left) and standard deviation (right) of superposition of glitch signals for different values of 0. Parameters: Af 
T = 0.01, <2>/(27r) = 1. 
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The total radiated energy per unit area per glitch is then 



dS 



dtdS 



IQtt^G 



dv 



T 


dt 


T 


\dh^,n 
dt J 



(A13) 



where the second equality is due to Parseval's theorem and 
T\\ denotes the Fourier transform. From Eq. (|A13|I . we iden- 
tify 
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as the time-averaged energy flux per unit frequency. The 
beam pattern, T.j. ' , which depe nds on the observer's po - 
sition through the angles 4> and t (jJaranowski et all 1 19981 ) . 
is ensconced in /ij^^. Averaging over line-of-sight orientation 
by integrating the right-hand side of E q. HA14|) with respect 
to cosijf) and i (jjaranowski et al.lll998l ). we can write 
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/;"#|r^'^--^2yifc,B2, ^^^ J^TT 



"g-i^) - i6^2Gp„r-^ -^ «. -^ ^. ' (A15) 

with i> = VT, /ts2 

We now calculate \T[dh^'i^ /dt]T[dh^'i^ /dt]\ for glitches 
arriving from an ensemble of neutron stars. This is a stochas- 
tic quantity, because h^^ depends on the random variables 
Ro and 0o- The mean and variance of Slgw are found in 
two stages. First, we calculate J^[9/iJ^/9i] for an individual 
glitch by applying the central limit theorem to add up the 
contributions from the AA'^v vortices involved. Second, we 
use properties of the x^ distribution to add up the signals 
from multiple glitches in multiple pulsars. 

According to the central limit theorem, the mean and 
variance of \T[dh^'k /dt\J^[dh^]^ /dt\\, for a glitch involving 
AA^v ~ sNvRs/Ar vortices, are ANvfii and AN^al respec- 
tively [fii and af are the Fourier transforms of Eq. p9p and 
Eq. (|20|) ]. The single- vortex moments fii and ai are calcu- 
lated by marginalising over the PDFs for initial vortex po- 
sition in Eq. (|17p and Eq. psp . the glitch size distribution 
g{s), and the source distance distribution A(D). 

The total \T[dhJiF /dt]T[dhJiF /dt]\ is the sum of con- 
tributions from Or simultaneous glitches. Its statistics obey 
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the non-central x^ distribution, with non-centrahty param- 
eter A-^2 = OrAA^vA*i/o"i, and mean and variance given by 
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QraJAN^ 1 -f 
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(A17) 



We note that, strictly, Eq. (|A16|) and (|A17|) apply only 
when AA'^v is constant. The correct expressions contain sums 
over different glitch sizes. However, comparison between 
Monte-Carlo simulations and Eq. (IA16|I and (IA17|I confirm 
that the approximation is reasonable. For simplicity, we also 
assume tj to be the same for all pulsars; the observed pop- 
ulation of glitching pulsars is dominated by young pulsars 
with 10 < a;/(27r) < 100 Hz. 

In the left panel of Fig. [XS] we graph hy, (t) from Monte- 
Carlo simulations for uj/(2ii) > 1 {top), uj/{2tt) = 1 (top 
centre) and a;/(27r) < 1 (bottom centre), representing vortex 
motions that are shorter, the same, and longer than one ro- 
tation period respectively. In the right panel, we graph the 
corresponding ilgw(i') (grey curves), overplotted with ana- 
lytic predictions of the mean (asterisks) and standard devi- 
ation (error bars). The error bars are too small to be seen. 
The bottom panels display results for a;/(27r) = 0, where 
radial vortex motion is the only contributor to the gravita- 
tional wave signal. 

A fair comparison between simulation and theory is 
achieved by restricting the glitch size to range between 
5 X 10~^ and 1 x 10~®, to ensure that the Monte-Carlo 
algorithm comprehensively samples g(s). In all four rows 
in Fig. IA3I the analytic mean is an excellent match to the 
Monte-Carlo results, but the analytic standard deviation un- 
derestimates the simulated value by a factor of ~ 100. For 
(l}/{2n) > 1, f2gw(i') peaks aX v = uj/(2-n) and is oscillatory, 
with nodes at integer multiples of u, except at i/max = Oj. For 
uj 1(211) < 1, r2gw(£') peaks at z>max ~ 0.375. The latter result 
holds for all uj/(2tt) <^ 1, because the canonical parameters 
listed in Table [T] describe a glitch in which vortex motion 
is dominated by azimuthal rather than radial motion. The 
peak at Anax ~ ijj/(2tt) is roughly five times higher than the 
next highest peak. 

We estimate rjgw as 
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where we use s+ = Af(/s//c) and s_ = Af/Nv These re- 
sults should be compared to estima tes of the backg round 
from sources of cosmological origin l|Maggiorel |2000| ) , such 



as cosm ic strings [Og 



lo- 



in the frequency ra nge 1.5- 



.w ~ 10" 



(|TurneiJ 
< 10-12 



2.5 kH z l|Ferrari et al.lll999l ')]. inflation [n 
[l997|)] and core-collaps e supernovae [10' 
at 100 Hz (|Buonanno et al. 2005)]. 

Although the background looks relatively weak at first 
blush, it could be much larger than the canonical numbers in 
Eq. HA19|I if (1) r is smaller than 10"^ s, which is eminently 
possible, and (2) extragalactic pulsars are included, since one 



has Q ex (-D)"^ as discussed previously. These possibilities will 
be explored in future work. 



A4 Signal to noise 

The optimised sign al-to-noise ratio (S/N) f or an integration 
time T is given by (|Allen fc Romandll999l ) 



S_ 

N 



•dH^T 



50714 



dv 



'i\v)nl.„(v) 



'''Pl(l^)P2(u) 



1/2 



(A20) 



where Pi(y) and P2(i') are the power spectral noise densi- 
ties for the two detectors, and 7 is the normalised overlap 
reduction function, which characterises the loss of sensitivity 
due to the separation and relative orientation of the detec- 
tors. The minimum observable figw of two detectors with the 
sensitivity of the Einstein Telescope (with a signal-to-noise 



ratio of 2.56) is Q,n 



1.13 X 10 for an integration time 



of one year, a false alarm rate of a — 10%, and a detection 
rate of 90% l|Marassi et al.ll2010l ). For a;/(27r) = 200 Hz, 
and assuming a Milky Way distance distribution given by 
Eq. (|A10|) . with an aggregate glitch rate of O = 10^ s-^, 
and glitches lasting r — 0.01 s, we obtain S/N — lO"'' for 
Milky Way glitches (with 7 = 1 for simplicity) . Very crudely, 
assuming there are ~ 10^" Milky- Way type galaxies, evenly 
distributed throughout the universe, S/N oc f2gw oc (D) be- 
comes ~ 1. 
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Figure A3. Superposition of gravitational wave signals in the cross polarisation hx (t) (in units of Ko/r"^; left column) and gravitational 
wave energy density parameter, Q.g^{v) {right column), from a population of glitching pulsars with mean aggregate glitch rate = 10, 
for cj/(27r) = 2.0, 1.0, 0.1 and (top to bottom respectively). The grey curves in the right column are different realisations of the 
underlying glitch size and distance distributions. The asterisks correspond to the theoretical calculation described in Sec. |X] Minimum 
and maximum glitch sizes arc 5 X 10~^ and 10~® respectively. Unless otherwise stated, the pulsar parameters arc given in Tabic, [l] 
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